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The crystal structure of boron is unique among chemical elements, highly complex, 
and imperfectly known. Experimentalists report the /3-rhombohedral (black) form 
is stable over all temperatures from absolute zero to melting. However, early calcu- 
lations found its energy to be greater than the energy of the a-rhombohedral (red) 
form, implying (5 cannot be stable at low temperatures. Furthermore, (5 exhibits 
partially occupied sites, seemingly in conflict with the thermodynamic requirement 
that entropy vanish at low temperature. Using electronic density functional theory 
methods and an extensive search of the configuration space we find a unique, energy 
minimizing pattern of occupied and vacant sites that can be stable at low temper- 
atures but that breaks the /3-rhombohedral symmetry. Even lower energies occur 
within larger unit cells. Alternative configurations lie nearby in energy, allowing the 
entropy of partial occupancy to stabilize the /3-rhombohedral structure through a 
phase transition at moderate temperature. 

PACS numbers: 61.66.Bi 

I. INTRODUCTION 

Elemental boron is important for its light weight, high strength, high melting point and 
semiconducting properties. It is also intrinsically interesting owing to its complex struc- 
tures which are characterized by their arrangements of icosahedral clusters. Boron is almost 
unique among elements in lacking a well understood and universally agreed upon low tem- 
perature structure. Knowledge of the precise structure is required for understanding its 



2 



remarkable electronic and mechanical properties, and also for understanding the stability of 
technologically important boron compounds relative to their constituent elements. 

The first^ reported crystallographic refinements of elemental boron were the a- 



^3, |4j. Although the a-rhombohedral form 



tetragonal and a-rhombohedral forms 

was initially believed stable at low temperature 3| , the /3-rhombohedral form, was later dis- 
covered j|| and proposed as the true low temperature state. Still later, the /3-tetragonal form 
was discovered and its crystal structure refined 6|. The a and /3-rhombohedral structures 
are illustrated in Figs. [1] and [2] and are described in greater detail in Section II Al 

The debate over the proper stable form of boron continues today, with some researchers 



suggesting a-rhombohedral is the true low temperature state 
energy of a lies below the reported calculated energies of (5 [7 



7]. In particular, the calculated 



10j. Some researchers 



propose a finite temperature phase transition Uj, includingone proposal that vibrational 



entropy drive, an „ t o , t_ at finite t e m perat„re ft Q U aat Um — zero 
point vibrational ener gy h as been proposed as a mechanism to stabilize the /3-rhombohedral 
form at absolute zero 10] . 

We note that the /3-rhombohedral form exhibits intrinsic disorder in the form of partially 
occupied sites. Partial occupancy occurs in crystallographic refinement when the site is 
occupied, at a given instant, in some copies of the unit cell but not in others. Within a single 
unit cell the site may be occupied at some times but not at others. Since thermodynamics 
requires that entropy vanish in the limit of low temperature, and partial occupancy implies 
finite entropy, the partially occupied /3-rhombohedral form is not a plausible low temperature 
structure. Correlations among the partially occupied sites of (3 must favor a unique pattern 
of occupancy that minimizes the energy. 

We carry out an extensive study of the configuration space that explores specific res- 
olutions of partial occupancy, assigning atoms or vacancies to specific sites. Our results 
indicate that a particular symmetry-broken form of (5 achieves energy lower than a and thus 
is the true low temperature structure. A symmetry restoring phase transition, driven by the 
high entropy of an ensemble of nearly degenerate configurations, should occur at a moderate 
temperature, explaining the experimental observation of /9-rhombohedral as the equilibrium 
phase at higher temperatures. 

The remainder of this introduction briefly describes the known structures of crystalline 
boron and our calculational methods. Results are presented in Section [III We find optimal 
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assignemnts of atoms to partially occupied sites, we validate these findings using solid state 
molecular dynamics simluations, and we explore possible superlatttice ordering. Finally we 
discuss our conclusions in Section IIHI 



A. Structural description 

Crystalline boron occurs primarily in tetragonal and rhombohedral forms. Since the 
tetragonal forms are known to be metastable, our analysis in this paper focuses on the 
rhombohedral forms. Both the a and /3-rhombohedral forms share a common space group, 
R3m (group # 166). Three rhombohedral primitive cells combine to form a larger nonprim- 
itive hexagonal unit cell with three times the number of atomic sites. The rhombohedral 
(111) axis becomes the hexagonal symmetry axis of the hexagonal unit cell. For simplicity 
our analysis is based on the rhombohedral primitive cell, though we do at the end consider 
supercells of the primtive cell including the full hexagonal unit cell. 

Structure types are denoted by their Pearson symbols. For example, Pearson type tP50 
(a-tetragonal) is primitive tetragonal with 50 atomic sites per cell. Pearson tP196 {(3- 
tetragonal) is also primitive tetragonal but with 196 atomic sites per unit cell. Because 
of partial occupancy tP196 actually has fewer than 196 atoms distributed among its 196 
sites. Pearson hR12 (a-rhombohedral) and hR105 (/^-rhombohedral) are both rhombohedral 
primitive cells of 12 and 105 sites repsectively. Three rhombohedral primitive cells may be 
combined to form a single hexagonal cell. 

All crystalline boron structures share a common structural motif, the 12-atom icosahedral 
cluster. Various allotropes differ in the spatial arrangement of the icosahedra and in the 
presence of interstitial atoms. The structural complexity of boron is due to the "electron 



deficiency" of the B12 icosahedral cluster 12J, [13] which frustrates the distribution of elec 



trons among available bonds. Presumably the partial site occupancy serves to relieve this 
frustration [14]. 

Among boron allotropes, a-rhombohedral boron has the simplest structure (see Fig. [1]), 
with a B 12 icosahedral cluster placed at each vertex of the rhombohedral cell. The structure 
is defined by just two independent Wyckoff positions [l5]. The icosahedral clusters are 
nearly regular, with bond lengths in the range of 1.75-1.81 A. Clusters are joined along six 
of their twelve 5-fold axes that point radially outwards through vertices of the icosahedron. 
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These intericosahedral bonds have length 1.67 A. Also visible in the figure are bonds of 
length 2.01 A running parallel to the icosahedral 2-fold axes. This bond length is almost 
completely absent in all other allotropes, so we believe it is energetically unfavorable. 

The /3-rhombohedral structure (see Fig. [2]) differs from a in that the icosahedral cluster 
at the origin joins to twelve other icosahedra along each of its twelve five-fold axes. In 
contrast to a, where all icosahedra share a common orientation, in (3 the twelve surrounding 
icosahedra are each rotated by 36° relative to the central icosahedron around their common 
axes. Consequently, the lattice constant of j3 is approximately double that of a. As shown, 
the icosahedra shaded in blue are at positions equivalent to those in the a structure, while 
those shaded in red are new. 

Another new feature in (3 is the presence of partially occupied sites. Although the orig- 
inal hR105 structure model 0, [if]] contained 15 Wyckoff positions, each fully occupied, an 



intermediate model 



17|,Q (Pearson type hRlll) assigned the B13 position 73% average 



occupancy and introduced a new B16 position at 25% occupancy (the numbe ring scheme we 
use is common to all the cited authors). In the Pearson type hR141 model 19( illustrated 
in Fig. El the B13 position is listed at 74.5% average occupancy, the B16 position at 27.2%, 
and additional positions B17, 18, 19 and 20 are given occupancies of 8.5%, 6.6%, 6.8% and 
3.7% respectively. Other than the B20 position, all partially occupied positions lie in mirror 
planes of the structure that contain the rhombohedral (111) axis. They appear as collinear 
sets in the projections shown in Figs. [2]and[3j 



B. Methods 



Our basic calculational methods follow Ref. [8j. The calculations use the electronic den- 
sity functional program VASP [20|] version 4.6.28. VASP uses a plane-wave approach that 
relies on periodic boundary conditions and is well suited to the study of periodic crys- 
tal structures. Electron-ion interactions are represented using Projector Augmented Wave 



(PAW) potentials |2ll . |22| which are an all-electron generalization of pseudopotentials. For 
boron the Is electrons are treated within the ionic core while the 2s and 2p electrons are 
assigned to the valence band. All structures considered were electrically neutral. 

The exchange-correlation functional is taken as the PW91 Generalized Gradient Ap- 



proximation (GGA) 23]. Previously the GGA has been shown superior to LDA (the Local 
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Density Approximation) for studies of Boron clusters and compounds [24], |25| , with accuracy 
nearly that of Hartree-Fock calculations at the self-consistent-field (SCF) level 26]. 

All structures are fully relaxed, both in atomic coordinates and lattice parameters, subject 
to the preservation of initial symmetry, using a conjugate gradient algorithm. Owing to 
coupling of the basis set to the volume we perform consecutive calculations (i.e. stop and 
restart) to fully relax the structure. Reported energies are obtained from a final static 
calculation. 

We test the dependence of energy differences among the a, /3-rhombohedral and optimized 
symmetry-broken (3 (see section IIIBj) structures on computational parameters. The results 
are summarized in table [TT1 To maintain consistent cell sizes we use a 2x2x2 supercell 
of the a phase, whose lattice parameters become similar to those of f3. The plane wave 
energy cutoff E cut is 319 eV by default, and we also test 415 eV. "Precision" settings of 
"Medium" and "Accurate" set the density of FFT grids which control wrap-around errors. 
We systematically increase the Monkhorst-Pack Appoint mesh starting from lxlxl (i.e. the 
T-point) until sufficient convergence is achieved. Based on the data presented we claim 
energy differences are converged to within 6% (roughly two significant figures) provided we 
use a fc-point mesh of 3x3x3 and "Accurate" precision. The default plane wave energy cutoff 
of 319 eV suffices. 

In Table II III we test the dependence on choice of potential and exchange-correlation 



functional. The local density approximation [27 



>DA) is expected to be the least accurate 
28] (USPP) treats only the valence electrons 



method. The ultrasoft pseudopotential method 
explicitly. The "hard" potential has a very small core radius and correspondingly high 
E cut =700eV . Results with this potential are expected to be comparable in accuracy to all- 
electron FLAPW calculations and to Gaussian with large basis sets, according to VASP 
documentation |29| . 

Because our final energy difference E opt — E a is not much larger than the variation 
among the GGA calculations it would be desirable to repeat this calculation using more 
refined quantum-chemical methods. However, all methods agree that our optimized struc- 
ture achieves a lower energy than a, with the exception of LDA which is expected to be the 
least accurate. Even if a more accurate calculation were to find that all variants of (3 were 
higher in energy than a, that would not alter our central conclusion that the fully symmetric 
/3-rhombohedral structure is a high temperature phase stabilized by occupancy fluctuations. 
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II. RESULTS 

Total energy calculations depend on precise knowledge of atomic positions, with partially 
occupied sites resolved into a specific pattern of occupied or vacant positions. Likewise 
partial site occupancy is thermodynamically forbidden in the T = OK limit. Hence we ex- 
plore the ensemble of likely instantaneous configurations, seeking both the unique optimal 
arrangement of atoms among partially occupied sites as well as an estimate of the entropy as- 
sociated with nearly optimal configurations. We carry out the study initially within a single 
rhombohedral primitive cell, repeated infinitely owing to the periodic boundary conditiions. 
Later we study superlattice ordering within supercells. 

TableUlists our main results [30]. Energies are given relative to the a-rhombohedral form, 
which we take as a reference because it contains no partial occupancy and was previously the 
lowest energy structure that was known. As expected, the tetragonal structures, known to 
be metastable, exhibit relatively high energies. The fully occupied hR105 ,5-rhombohedral 
structure is higher in energy than the a-rhombohedral hR12 structure, as previously noted 17, 

nnn 

|8|, 121, HOj. If (3 is to be stable at low temperatures it must involve the placement of atoms 
among the partially occupied sites of hRlll or hR141. Some particular pattern of occupied 
and vacant sites must minimize the total energy while breaking the rhombohedral symmetry. 

The following tables present energy data supporting specific conclusions on the optimiza- 
tion of B. hRlll and B.hR141. Figures and atomic coordinates for each named structure, 
and many other structures not listed here, can be viewed at our web site {scj (see the special 
"published" area). Data is shown for symmetry- inequivalent structures. All energies are 
given for a 3x3x3 Monkhorst-Pack fc-point mesh, for accurate precision and the standard 
PAW potential with default plane wave energy cutoff of 319 eV as discussed in Section II Bl 

The notation lists only those sites among the partially occupied Wyckoff positions that 
are actually occupied. Positions B1-B12, B14 and B15 are always fully occupied and thus 
are not listed in our notation. The hRlll structure introduces partial occupancy at B13 
and the new B16 position. The hR141 structure introduces additional positions B17-B20. 
The notation abode f refers to specific sites within each Wyckoff position, as labeled in Fig[3l 
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A. Pearson hRlll 



We first explore the hRlll model [18j. The two partially occupied sites, B13 and B16 
are both of Wyckoff type 18h, meaning that 18/3=6 of these sites occur per rhombohedral 
primitive cell. The B13 sites form a pair of equilateral triangles surrounding the B15 site 
at the primitive cell body center (see light blue atoms in Fig. [3]). The B16 sites form a pair 
of equilateral triangles that lie immediately above faces of the icosahedra at primitive cell 
vertices (see pink atoms in Fig. [2]). 

A B13 atom can be swapped for a B16 atom in 4 symmetry- inequivalent ways, each of 
which lowers the energy by 12-13 meV/atom. This substitution is thus strongly preferred 



energetically, but the spatia 
weak, as suggested by Slack 



correlation between occupied and vacated site is relatively 



19|]. 



Since the number of atoms per cell is believed to be greater than 105, we also considered 
a single B13 vacancy and a pair of B16 atoms, again exhaustively testing all combinations. 
Low energy requires that one atom reside on each equilateral triangle of B16. Our optimal 
structure within the confines of hRlll is B13bcdefB16bd, at an energy of 0.15 meV/atom 
above hR12. The B13bcdefB16bf structure was the best found in a previous study by van 
Setten, et al [10 ]. 



B. Pearson hR141 



So far the occupancy of B13 is larger than reported experimentally, the energy remains 
above the energy of a, and the total number of atoms per primitive cell remains below the 
experimentally reported value of 106.7. Evidently we should try removing another atom 
from B13 and placing that atom (and more) at other locations. To place these atoms we 
utilize the additional partially occupied positions of hR141, namely B17-20. There are of 
order 10 6 distinct arrangements of atoms within a single primitive cell that are consistent 
with experimentally observed occupancies. Since this is far too many structures to explore 
exhaustively, we build upon our prior results and consider other likely correlations in order 
to focus our search. 

n 

Following Slack [19] we note the 1.57 A bond between B13 and B17 sites is slightly 
too short for simultaneous occupancy. If we choose to occupy B17a, this model suggests 
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that B13d and B13e should be vacant. However, our calculations show this arrangement 
is not stable. One of the vacancies moves to the B13a site leaving just a single vacancy 
on either B13d or e but not both. Then the B17a atom relaxes to accommodate the bond 
to the remaining nearby B13 atom, which ends up at length 1.86 A. It is noteworthy that 
the B17 site reports an anomalously large thermal Debye-Waller factor, indicating large 
displacements from the refined position. We suggest that this position should be split into 
distinct sites whose occupancy is correlated with the nearby B13. In fact, Slack utilized an 
alternative split site "B17d" to refine one of his samples, and one of these B17d sites lies 
within 0.1 A of our relaxed B17 position. 

The 1.62 A bond length between sites B17 and B18 is within the favorable range, and 
their reported occupancies are similar, so we presume their occupation is correlated. Our cal- 
culated energies confirm this, with a reduction in energy of 1.1 meV/atom upon introduction 
of an atom at the B18 site adjacent to an occupied B17 site. After optimizing the placement 
of B16 atoms this leads to our optimal primitive cell structure, B13bcefB16bdB17aB18a, 
which contains 107 atoms and achieves an energy 0.86 meV/atom below the energy of a. 
This model also achieves a close match of relaxed positions to experimentally observed po- 
sitions, as is evident from the AR values in Table [B 

The B19 and B20 sites so far remain unused in our study. B19 sites are only 0.74 A away 
from B18 sites, so these may never be simultaneously occupied. Both B19 and B20 sites lie 
close to the centers of fully occupied hexagonal rings, resulting in small patches of triangular 
lattice that are atypical of crystalline boron structures (although they may be stable in small 
boron clusters and nanotubes 31|). We found no significantly low energy structures utilizing 



B19 or B20 sites. Van Setten et al. 



lOj report a low energy for B19bcdefB16eB19a, but we 



find this is not favorable. Rather the B19a atom relaxes to the nearby B17a position. 



C. Molecular dynamics 

Because we cannot systematically evaluate all configurations within the hR141 primitive 
cell, we checked our result using molecular dynamics. We used the VASP-calculated forces 
to perform molecular dynamics simulations of a single primitive of (3. To achieve atomic 
diffusion we employed a high temperature T=2000K (melting is around T=2365K). We 
estimated the lattice parameters as a= 11.047 A andc=24.155 A based on an experimental 
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report of thermal expansion 32]. 

We ran samples of 106, 107 and 108 atoms for a duration of 16ps each, using lfs time 
steps. To search for optimal configurations we drew instantaneous configurations from the 
molecular dynamics run every 2ps and quenched them. To carry out the quench we rescaled 
the lattice parameters to their crystallographic values and performed molecular dynamics 
with a linearly decreasing temperature ramp that reached T=300K after 3ps. We then 
performed conjugate gradient relaxation to reach OK. This procedure was able, on occasion, 
to achieve the optimal structures we reported, but usually resulted in higher energies. 

Fig. [5] shows a density plot of atomic positions for the 107 atom run that started with 
our optimal structure as an initial condition. Densities have been averaged to impose rhom- 
bohedral symmetry. Sharing of atoms between B13 and B17 positions is clearly visible, as 
is partial occupation of B18 and B16 positions. 



D. Supercell studies 

Although we optimized the assignment of partial occupancy within a single hR141 prim- 
itive cell, there is a possibility of a lower energy structure within either the hexagonal unit 
cell (Pearson type hP423) or some other supercell of the primitive cell. 

The hexagonal unit cell (Pearson type hP423, dimensions a=10.9 A, c=23.8 A) contains 
three copies of the rhombohedral primitive cell stacked in the direction parallel to the rhom- 
bohedral (111) axis. We consider the case where each primitive cell is identically decorated, 
where one is rotated by 120° relative to the other two, and where each is rotated by 120° 
relative to its neighbors. In this last case a 3i screw axis is introduced. In Table PVTl energies 
are given for a 3x3x1 A;-point mesh with medium, maintaining a uniform reciprocal space 
density comparable to a 3x3x3 mesh for a single primitive cell. It appears energetically 
preferable to maintain identical orientations of all vertically stacked cells. 

However, within a 2x1x1 supercell of the 107 atom primitive cell, which places indepen- 
dent primitive cells adjacent to each other (see Fig. HJ), we did find superstructures which 
lowered the energy significantly. To describe these we focus on a motif near the center of our 
optimal primitive cell structure. One of the two B13 vacancies is collinear with the occupied 
B17 and B18 sites (e.g. B13aB17aB18a). The second B13 vacancy is adjacent to the first 
(e.g. B13d). This motif of collinear and adjacent sites can occur in six rotated variants and 
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six more that are reflected versions. 

We examined all 22 symmetry- inequivalent arrangements of this motif within the 214 
atom supercell. For each arrangement we optimized the placement of the B16 atoms. The 
relative arrangements are coupled mainly by the placement of B16 atoms, and all energies 
lay within 2 meV/atom of each other. Energies are given in Table lyTTI for a 1x3x3 fc-point 
mesh and precision set to "Accurate". Several supercell structures yielded energies lower 
than the optimal 107 atom structure, indicating a preference for superlattice ordering at 
low temperatures, with adjacent cells of (3 resolving their partial occupancy with differing 
orientations of a common motif. 

The large number of nearly degenerate configurations suggests the possibility of a phase 
transition from the symmetry-broken low energy structure to a state that restores the (3- 
rhombohedral symmetry by sampling the full ensemble of motif orientations. To judge the 
chance of such a phase transition, we evaluate the partition function of our 2x1x1 supercell 

Z = J2^ a e- Ea/kBT (1) 

a 

where a runs over all 22 symmetry-independent configurations, Q a is the multiplicity of the 
configuration and E a is the relaxed energy, optimized over placements of B16 atoms (as in 
Table IVHl but multiplied by 214 for the number of atoms per supercell). 

Thermodynamic derivatives of Z yield the internal energy U, the entropy S and the 
heat capacity C. The heat capacity exhibits a strong peak around T = 300K. Fig [6] plots 
thermodynamic data resulting from this model. The strong peak around T=300K represents 
the unlocking of the relative orientations of our collinear motif (see above). The small 
heat capacity peak around T=50K represents the unlocking of the second B13 vacancies 
(see above) while holding the collinear motif fixed. The entropy available from unlocking 
dominates the energy cost, substantially lowering the free energy. Based on this it seems 
likely that in an infinite system there should be a symmetry-restoring phase transition at 
a moderate temperature, driven by the entropy of suboptimal occupation patterns of the 
partially occupied sites. This transition may be difficult to observe experimentally because 
atomic diffusion is slow at this temperature. 
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III. DISCUSSION 

The symmetry space group of /3-rhombohedral boron is R3m (group # 166). To pre- 
serve rhombohedral symmetry in the primitive cell, every Wyckoff position must be fully 
occupied or fully vacant. If the favored structure occupies only a subset of the sites in one 
or more Wyckoff position, the symmetry is necessarily lower. Our optimal vacancy-ordered 
structure within the primitive cell yields space group PI (group #1) corresponding to no 
point symmetry whatsoever. 



According to Landau theory [33( changes of symmetry occur through thermodynamic 
phase transitions, implying that the /5-rhombohedral structure is not the low temperature 
stable phase. Experimental observation of (3 as stable indicates either that vacancies are 
frozen in a disordered arrangement or that (3 is stabilized by a symmetry-restoring phase 
transition at some low temperature. This conclusion holds regardless of whether we have 
found the true optimal structure, since all energetically plausible variants of (3 lack full 
symmetry. 

The transition temperature is only a crude estimate because of several approximations. 
We consider only the coupling of primitive cells. We consider only a small fraction of all 
configurations within these two cells. We neglect atomic vibrations which contribute strongly 



to the free energy 0, 



lo| and must be included in any attempt to estimate an accurate 
transition temperature. However, while the a and (3 phases differ strongly in vibrational 
free energy, the ensemble of nearly optimal symmetry-broken f3 variants most likely remains 
nearly degenerate in vibrational free energy because their local environments are so similar. 

Because our calculated energy differences are at the borderline of reliability of density 
functional methods it would be desirable to study these energies using more sophisticated 
quantum chemical methods. Unfortunately at present these methods cannot be reliably 
applied to periodically repeated crystal structures of the necessary complexity. Although 
details of the optimal structures may change, our primary conclusions are unlikely to be 
altered by higher accuracy. The true ground state will be either symmetry-broken (3 or 
perhaps a. In either case a phase transition is necessary to recover the experimentally 
observed equilibrium state above some moderate temperature. 

Even within the density functional theory limitations there remains some uncertainty 
concerning the optimal structure. Within a single primitive cell we have only explored a 
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fraction of all configurations and some other unexamined structure might turn out to be 
prefered. Among supercells we have only addressed the hexagonal unit cell and the 2x1x1 
supercell. It is quite possible that a larger supercell will be even more favorable. In fact 
it could be that the true ground state even restores rhombohedral symmetry, though in a 
supercell so that it is no longer the f3 structure. 

In conclusion, we show that suitably resolving correlations among partially occupied 
sites of (3 yields an optimal structure whose energy is lower than a. This structure breaks 
the symmetry of (5 because it uses only subsets of the fully symmetric Wyckoff positions. 
Further, we demonstrate the likelihood of superlattice ordering. We propose that the full 
symmetry of (3 is restored at moderate temperatures through a phase transition driven by 
the entropy of partial site occupation. 
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Name 


Pearson 


Atoms 


V 


E — E a 


AR 


Comments 


a-T 


tP50 


50 


7.67 


91.91 


0.057 


full occupation 


(3-T 


tP196 


192 


7.53 


15.13 


0.056 


optimized occupation 


(3-R 


hR105 


105 


7.72 


25.87 


0.023 


full occupation 


a-R 


hR12 


12 


7.18 


0.00 


0.002 


full occupation 


(3-R 


hRlll 


105 


7.69 


13.02 


0.013 


B13bcdefB16a 


(3-R 


hRlll 


106 


7.64 


0.15 


0.012 


B13bcdefB16bd 


(3-R 


hR141 


107 


7.57 


-0.86 


0.005 


B13bcefB16bdB17aB18a 


(3-R 


hR141 


108 


7.56 


1.43 


0.015 


B13bcefB16acdB17aB18a 


2x1x1 


aP282 


214 


7.57 


-1.75 




optimized supercell 



TABLE I: Structural data including atoms per primitive cell, volume (A 3 /atom) and energy relative 
to a-R (meV/atom). AR (A) measures the deviation of the symmetry-averaged relaxed positions 
from the crystallographically reported positions, averaged over the fully occupied Wyckoff positions. 
Comments list occupied sites using notation in Fig. [3l From top to bottom: tetragonal, fully 
occupied rhombohedral, rhombohedral symmetry-broken structures and our optimized supercell 
structure. 



AE 




Ep- 


- E a 






E op t 


— E a 




E C ut 


319 eV 


415 eV 


319 eV 


415 eV 


Precision 


Med 


Acc 


Med 


Acc 


Med 


Acc 


Med 


Acc 


fc=lxlxl 


-7.76 


-7.54 


-7.98 


-7.97 


-17.27 


-17.08 


-17.04 


-17.15 


Zc=2x2x2 


23.58 


23.78 


23.85 


23.77 


-1.71 


-1.47 


-1.47 


-1.68 


A: =3x3x3 


25.71 


25.88 


26.09 


25.99 


-1.05 


-0.86 


-0.65 


-0.85 


&=4x4x4 


25.29 


25.39 


25.65 


25.49 


-1.04 


-0.88 


-0.64 


-0.87 



TABLE II: Test of convergence as function of energy cutoff, VASP precision setting and /c-point 
mesh. Energies of hR105 (Ep) and our optimal B13bcefB17aB18a structure (E op t) are compared 
with hR12 (E a ). 
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AE 


E[3 - 


- E a 




- E a 


Precision 


Medium 


Accurate 


Medium 


Accurate 


LDA 


47.83 


47.86 


15.48 


15.46 


USPP 


27.82 


27.56 


-0.25 


-0.45 


HARD 


36.06 


25.94 


-0.59 


-0.75 



TABLE III: Test of alternate potentials. All calculations are done with a 3x3x3 mesh at the default 
energy cutoff. USPP and HARD use the PW91 GGA. LDA and HARD use PAW potentials. 



B13 sites 


B16 sites 


Atoms 


E — E a 


a b c d e f 






105 


25.84 


a b c d e f 


a 




106 


12.25 


a b c d e f 


b 


d 


107 


12.50 


b c d e f 






104 


27.38 


b c d e f 


a 




105 


13.02 


b c d e f 


a b 




106 


13.81 


b c d e f 


b 




106 


13.10 


b c d e f 


b 


d 


106 


0.15 


b c d e f 


b 


e 


106 


5.21 


b c d e f 


b 


f 


106 


0.87 



TABLE IV: Table of selected data for hRlll. (top) Full B13 occupation; (middle) Single B13 
vacancy with single or double occcupation of one B16 triangle; (bottom) Single occupancy of both 
B16 triangles. Energy units are meV/atom. 
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B13 sites 


B16 sites 


B17 


B18 


B19 


B20 


Atoms 


E - E a 


a b c d e f 


b d 


a 








108 


6.62 


b c d e f 


b d 


a 








107 


20.29 


be e f 


b d 


a 








106 


0.15 


be e f 


a d 


a 


a 






107 


-0.07 


be e f 


b d 


a 


a 






107 


-0.86 


be e f 


c d 


a 


a 






107 


3.48 


be e f 


b d 






a 




106 


12.96 


be e f 


b d 


a 


a 




h 


108 


0.844 



TABLE V: Table of selected data for hR141. (top) B13 vacancies; (middle) B16 occupation; 
(bottom) Structures occupying B19 and B20 sites. Energy units are meV/atom. 



B13(l)(2)(3) 


B16(l)(2)(3) 


B17B18(1)(2)(3) 


Atoms 


AE 


(beef) (beef) (beef) 


(bd)(bd)(bd) 


(aa)(aa)(aa) 


321 





(beef) (acde) (beef) 


(bd)(cf)(bd) 


(aa)(bb)(aa) 


321 


0.09 


(beef ) (acde) (abdf ) 


(bd)(cf)(ae) 


(aa)(bb)(cc) 


321 


0.31 



TABLE VI: Hexagonal unit cell energies. Site occupancy is given for each primitive cell, (1), 
(2) and (3). AE is relative to a supercell of the optimal 107 atom hR141 structure, in units of 
me V/ atom. 



B13(l)(2) 


B16(l)(2) 


B17B18(1)(2) 


n 


AE 


(beef) (acde) 


(ad)(bd) 


(aa)(bb) 


8 


-0.89 


(bcef)(acef) 


(ad)(bd) 


(aa)(bb) 


4 


-0.75 


(bcef)(bcdf) 


(ce)(ad) 


(aa)(ee) 


8 


-0.24 


(acde)(bcdf) 


(bd)(ad) 


(bb)(ee) 


8 


-0.11 


(beef) (beef) 


(bd)(bd) 


(aa) (aa) 


12 


0.00 



TABLE VII: Selected 2x1x1 supercell energies, including multiplicities ft. AE is relative to a 
doubling of the optimal 107 atom hR141 structure, in units of meV/atom. 




FIG. 1: Structure of a-rhombohedral boron viewed along the rhombohedral (112) axis. Bond color 
scheme: 1.67 A in purple; 1.75-1.81 A in orange; 2.01 A in red. 
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FIG. 2: Structure of /3-rhombohedral boron viewed along rhombohedral (111) axis. Bond color 
scheme: 1.63-1.73 A in purple; 1.73-1.92 A in orange. Partially occupied sites (see text) are shown 
in color: B13 (74.5% average occupancy) cyan; B16 (27.2%) pink; B17 (8.5%) yellow; B18 ( 6.6%) 
indigo; B19 (6.8%) blue; B20 (3.7%) orange. 




FIG. 3: Structure of /3-rhombohedral boron in vicinity of cell body center (B15, shown in red) 
viewed along the rhombohedral (111) axis. Color coding as in Fig. [2J Site labels correspond to 
notation in table HI Size of atoms indicates vertical position, with small on top. Length scale is in 
Angstroms. The dashed lines contains a single hexagonal unit cell. 
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FIG. 4: Optimal site occupation of a /3-rhombohedral 2x1x1 supercell viewed along the rhombo- 
hedral (111) axis. Color coding as in Fig. [3l Small empty circles locate vacant sites. 
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FIG. 5: Density plot of simulated atomic positions. Grayscale indicates low frequency positions, 
colors indicate medium (blue) - high (red) frequency positions. Crystallographically determined 
atomic positions are superimposed (and suitably scaled) for comparison with Fig. [3J (top) Run 
started in optimal configuration, (bottom) Same data with rhombohedral symmetry imposed by 
averaging. 



22 




T (K) 

FIG. 6: Heat capacity (black), energy (red) and entropy (green) of supercell model. Energy is 
relative to the optimal supercell structure. 



